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Quantile regression provides a framework for modeling statistical 
quantities of interest other than the conditional mean. The regression 
methodology is well developed for linear models, but less so for non- 
parametric models. We consider conditional quantiles with varying 
coefHcients and propose a methodology for their estimation and as- 
sessment using polynomial splines. The proposed estimators are easy 
to compute via standard quantile regression algorithms and a step- 
wise knot selection algorithm. The proposed Rao-score-type test that 
assesses the model against a linear model is also easy to implement. 
We provide asymptotic results on the convergence of the estimators 
and the null distribution of the test statistic. Empirical results are 
also provided, including an application of the methodology to forced 
expiratory volume (FEV) data. 

1. Introduction. Quantile regression iias appeared as an alternative to 
least squares in a wide range of applications. When the center of the con- 
ditional distribution of a response variable Y , given a covariate vector X , 
is under investigation, median regression provides a consistent estimator of 
the conditional median without assuming a specific form for the conditional 
distribution. When other conditional quantities, for example, the lower or 
upper tail of the conditional distribution, are of interest, quantile regres- 
sion provides a way to directly estimate the interesting quantities without 
assuming that such quantiles are related to X in the same fashion as the 
conditional mean. A good example can be found in a study of consumer 
demand for electricity, where heavy users responded much more drastically 
to weather and time variation than average users [17]. 

Quantile regression has been implemented with nonparametric methods 
to overcome the limitation of a linear model. Local polynomial, smoothing 
spline and B-spline smoothing methods were considered in [4, 35], [13, 24] 
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and [14], respectively. He and Shi [15] used bivariate tensor product B- 
splines with a partly linear model. Wei and He [34] considered B-splines 
with a partly linear model for longitudinal data to model growth curves. 
De Gooijer and Zerom [6] and Doksum and Koo [7] considered local polyno- 
mials and polynomial splines with additive models. Nevertheless, compared 
to its development with a linear model, the nonparametric implementation 
of quantile regression has been limited. 

In this paper, we consider a varying-coefficient model for the conditional 
quantiles: for a random vector (T, X, Y) G [0, 1] x ]R^'+^ x M, we suppose that 

(1.1) qr{t, X) = Po{t)xo + ^l{t)xi + ■■■ + (3p{t)Xp, 

where Qrit, x) denotes the rth conditional quantile of Y, given that (T, X) = 
{t,x), and the are unknown smooth functions of t for j = 0, 1, . . . ,p. 

A random sample is modeled by 

(1.2) Yi = qr{Ti,Xi) + ei, l<i<n, 

where the Cj are independent random variables with rth quantile and are 
independent of the (Ti,Xi). 

Varying-coefficient models constitute an important class of nonparametric 
models: they have been widely applied to analyze conditional means due to 
their flexibility and interpretability (see [2], [9] and [19] for examples of the 
analysis of time series, survival and longitudinal data). However, a varying- 
coefhcient model for the conditional quantiles was considered only recently 
by [18] and [3] , where local polynomials were used for independent and time 
series data, respectively. In this paper, we provide a more comprehensive 
treatment of the model for independent data by proposing a polynomial- 
spline-based methodology for its estimation and assessment. 

The proposed methodology is quite general, in the sense that we do not 
require the unknown coefficients (3j{t) to be smooth curves of a common 
degree or the to be identically distributed; such assumptions facilitate 
the development of asymptotic theories covered in the later sections. The 
methodology is also simple and easy to implement: the proposed estimator 
is readily computed by standard quantile regression algorithms and the pro- 
posed stepwise knot selection algorithm. Different knots can be selected for 
different coefficients and the algorithm can accommodate varying degrees 
of smoothness without complicating the computation (see [11] and [19] for 
discussion of this issue with respect to local polynomials and smoothing 
splines). We show that the spline estimators and their derivatives attain the 
optimal rates of global convergence under appropriate conditions. 

The stated feature of accommodating different degrees of smoothness in 
the coefficient functions is shared by any estimation method employing poly- 
nomial splines, such as the one which appears in [21]. However, the one in 
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[21] requires its knots to be equispaced and selects the numbers of knots 
for the coefficients via cross-validation. In addition, it uses the least squares 
method for estimation since it concerns the conditional mean. 

For the assessment of the model against a linear model, a Rao-score-type 
test is proposed. The test uses polynomial spline smoothing to account for 
the structure of the varying-coefficient design space and is therefore poten- 
tially more powerful than a test with a nonspecific nonparametric alterna- 
tive. Implementation of the test is easy, since the test statistic has a simple 
asymptotic null distribution. 

Model assessment is essential in this context since a linear model will be 
preferred on the grounds of parsimony unless the underlying science war- 
rants use of the varying-coefficient model. While this issue is generally im- 
portant for nonparametric models, it has received little attention in quantile 
regression. Only recently, He and Zhu [16] and Horowitz and Spokoiny [20] 
considered general lack-of-fit tests for linear quantile regression. Both tests 
are consistent for any fixed alternative. However, they are not powerful in 
this context because their alternative space does not take into account the 
structure of the varying-coefficient model. Any lack-of-fit test with a non- 
specific alternative will have the same problem and is thus not considered 
here (see [12] for a survey of such general tests for the conditional mean). 

Likelihood-ratio-type tests were proposed for the conditional mean in a 
similar context [1, 2, 10]. Local polynomials were used to smooth the varying- 
coefficient alternative design space. In essence, these tests compared the 
residual sums of squares under the null and the alternative. We adapted 
and used the tests to evaluate the performance of the proposed test in a 
Monte Carlo study. The results are presented in Section 4.1. 

We present a motivating example involving forced expiratory volume 
(FEV) data collected on children aged three to nineteen [29]. A topic of 
interest is the effect of cigarette smoking on FEV after accounting for the 
effects of age, sex and height. One may speculate that the effects of the inde- 
pendent variables change with age. On the other hand, the lower conditional 
quantiles of FEV are of independent interest as they can provide a gauge 
for poor pulmonary function. Both research equations can be investigated 
by fitting model (1.2) and testing for age dependency in the effects of the 
independent variables. An important point is that we do not need to assume 
that the effects of the independent variables are same across the quantiles. 
Detailed results for this example will be given in Section 4.2. 

In this paper, asymptotic results will be proven under rather simplified 
conditions that require knots for the splines to be nonstochastic and quasi- 
uniformly placed. As with most other asymptotic analyses in the literature, 
the results obtained under such conditions can provide insight into the large- 
sample behavior of the proposed methodology in more general settings. The 
rest of the article is organized as follows. Sections 2 and 3 discuss model 
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estimation and assessment. Section 4 presents empirical results. Proofs are 
deferred to the Appendix. 

2. Estimation. 

2.1. Regression splines. The degree of smoothness of the true coefficient 
functions determines how well the functions can be approximated. We as- 
sume that the (3j (t) are functions of t with a common degree of smoothness 
r defined as follows: let TCr be the collection of all functions on [0, 1] for 
which the mth order derivative satisfies the Holder condition of order 7 
with r = m + 7. That is, for eachheHr, \h^"''> {s) - h^'^'> {t)\ < Wols-t^ for 
any < s, t <1 and a positive finite constant Wq. 

Condition 1. /3j{t) G Hr, j = 0, 1, . . . ,p, for some r > 1/2. 

If the have bounded dtli order derivatives on [0, 1] , then Condition 1 
holds with r = d. 

We shall use normalized B-splines of order m + 1 to approximate the f3j{t). 
We consider a sequence of positive integers {kn}, n > 1, and an extended 
partition of [0,1] by fc„ quasi-uniform knots. Following [30], we denote the 
associated B-spline basis functions hy Bj{t), j = 1, ... ,kn + m + l. We define 
TTfc^W = (-Bi(t),...,Sfc„+m+i(t))'^ and Ilk,,(t,x) = (3;o7rfc„(t)'^,a;i7rfc„(t)'^, . . . , 

The proposed polynomial spline estimator of qr{t,x) is given by 

p 

(2.1) (t, x) = Y, ej^u„ {t)xj = elu,^ (t, x) , 

j=0 

where pk„ = {p+ 1)(A:„ + m + 1) and = (^o^, • • • , OjV ^ solves the 
minimization problem 

(2.2) min ^ p^iY, - 0'^Uk„iT,,X,)) 
em^^n ^^^^ 

for Pr{s) = s{t — I{s < 0)). Accordingly, the polynomial spline estimator 
Pj{t) of Pj{t) is given by 9j7rk„{t) for each j. 

2.2. Asymptotic result. The following are sufficient conditions for the 
proposed polynomial spline estimator and its derivatives to converge at their 
best possible rates as the sample size goes to infinity. For simplicity, we adopt 
the following notation throughout the paper: for a vector v, denotes its 
jth element and \u\ denotes its Euclidean norm. We use a„ ~ bn to mean 
that there are constants < ^ < S < cxd such that A < fln/^n < B for all n. 
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Condition 2. The conditional distribution of T, given X = x, has a 
bounded density /t|jc : < ci < fT\x{A^) < C2 < oo uniformly in x and t for 
some positive constants ci and C2- 

Condition 3. X(o) = 1, E{Xf^j)\T) = and P{|X(j)| < M} = 1 for some 
M < oo, J = 1, . . . ,p. There exist two positive definite matrices Si and S2 
such that Si < Var(X|r) < S2 uniformly in T, where Var(X|T) denotes 
the conditional covariance matrix of X given T. 

Condition 4. The ej are i.i.d. and have a density function /g that is 
continuous at with < /e(0) < 00. 

With an appropriate choice of kn to balance bias and variance, the pro- 
posed polynomial spline estimator and its derivatives attain the optimal 
convergence rates established by [31], as shown in the following theorem: 

Theorem 1. Assume Conditions 1-4- Suppose fc„ ~ n^/^^^"^^) and r > 
1/2. Then for j = 0,1,..., p, 

(2.3) i )(r,) - pf\T,)f = 0,(n-2(^-'=)/(2'-+i)), fc = 0, 1, . . . ,m. 

i=l 

If piecewise linear splines are used and the (3j{t) have bounded second 
order derivatives, then Theorem 1 gives the rate of convergence as n~^/^. 

Remark 1. Considerable effort has been directed at relaxing the 'non- 
stochastic quasi-uniform knots assumption' in least squares regression. For 
example, Mao and Zhao [27] and Stone and Huang [33] investigated free- 
knot splines in which both the knot locations and spline coefficients were 
treated as unknown parameters. Extending the methodology and theory of 
free-knot splines to quantile regression is beyond the scope of this paper. 

2.3. Implementation. When the B-spline basis is given, computations 
can be performed using standard quantile regression algorithms. As for 
selecting the order and knots for the splines, we propose a simple semi- 
automatic stepwise algorithm: users determine the order and specify a set of 
potential knots and the knots are chosen automatically from the potential 
knot set. The algorithm adds or deletes knots iteratively using Rao and Wald 
statistics to avoid additional model fitting. It allows different knots to be 
chosen for different coefficients. The algorithm terminates when no further 
addition or deletion occurs or when the model at the previous iteration is 
the same as the one at the current iteration. Similar algorithms are found in 
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[26] and [32]. Upon termination, the algorithm finds the best fitting model 
using an adapted Schwarz-type Information Criterion (SIC), 

(2.4) SIC{M^)=log(^J2Prirj)^ +0.51og(n)K/n, 

where M* denotes the ith interim model, rj is the jth residual from fitting 
M*, pI^ is the number of variables in M* and n is the sample size. For 
algorithmic convenience, we use a truncated power basis. We refer to Chapter 
3.2 of [22] for more details. 

2.3.1. Order of splines to fit. We suggest using lower order splines, such 
as linear (m = 1) and quadratic splines {m = 2), for practical reasons. Since 
the effect of the splines on the model is multiplicative, higher order splines 
would induce complicated interactions and collinearity among the variables 
in the model. For example, even the simplest cubic splines (one-piece cubic 
polynomials) would induce interactions of the form Xjt, Xjt^ and Xjt^ . We 
prefer linear splines in particular since with them, the coefficient estimates 
are easier to interpret. Moreover, linear splines have an optimal property 
[24]. 

2.3.2. Potential knot set. Users provide either a potential knot set or 
the number of potential knots. In the latter case, the algorithm chooses 
as potential knots the user-provided number of equispaced knots. A helpful 
guideline of min(4n^/^, n/4, A^, 30) is found in [32] for the number of potential 
knots, where is the number of distinct data points. 

3. Model assessment. Model assessment takes the form of a hypothesis 
test examining whether all of the coefficients are constant. Under Condi- 
tion 1, this hypothesis test can be represented as 

(3.1) Ho:qr{t,x)£C vs. Hi:qrit,x) e g\C, 
where C and G are sets of functions defined as 

£ = {q{t, x) I q{t, x) = boXQ + • • • + bpXp, for some constants bj, j = 0, . . . ,p}, 

Q = {q{t,x) I q{t,x) = (3o{t)xo H \- (3p{t)xp, for (3j{t) G 7ir, J = 0, . ..,p}. 

We consider a transformation of Ilf:^{t,x): Ti.^{t,x) = Afc^nfc^(t,x) = (x"^, 
r2,fc„(t,x)"'")^ for some transformation matrix A^^. We denote the vector 

of coeflacients of Tk„{t,x) by i^r, = iii,kn^^2,kj'^ ^ ^^^^^ ^i,k„ and la.fc^ 
correspond to x and T2^kn (t^x), respectively. Then the minimization problem 
is 

n 

(3.2) min^ Pri^^ ^Ijk„m,X,)) 
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and (3.1) can be represented as Hq : $2,k„ = versus Hi : ^2,fc„ 0- We denote 
the estimates of ^^.^ obtained under Hq and Hi by ^/.^ = (^1^,0^)"'" and 

^fc„ = (IZfc„'U,fcJ^' respectively 

The proposed Rao-score- type test uses the score of (3.2) evaluated at the 
estimates of^^^ undergo. Let Sk„ = {kn/n)^/'^J2i=i^riYi-iJXi)r2,kniTi,Xi), 
where (/7r(e) = r/(e > 0) + {t — l)/(e < 0), the derivative of Prie)- Define 
a^ = E{ip^{ef) and Q„ = Er=irfc„(T„X,)rfc„(T„X,)T. Let Q„(,^), i,j = 

1,2, denote the ijth block of Q„ such that Q„(i2) = X^iLi ^jr2,fc„ (T,, Xj)"^, 
for example. The test statistic is 

(3.3) ^fc„ = (l/a)(fe„Q(22)/n)i/2,^^, 

where Q^^^^ = (Q„(22) - Qn(21)Qn{ll)Qn(12))"^- 

The test statistic has a simple asymptotic null distribution. Suppose that 
kn = kQ, a positive constant. Under Conditions 1-4, for any qr{t,x) £ C, 
l^fcoP — X(p+i)(fcp+i) as n — > oo [23]. For increasing kn, we have the fol- 
lowing theorem: 



Theorem 2. A ssuTfic CoTiditioTis Suppose lini7^_>(2>o kln'-^ = for 
some < (5 < 1. Then for any qr(t,x) G C, 

— > A'(0,1) as pk„^ CO. 



V2(Pfc„-p-l) 

If kn is bounded, then one uses the chi-square distribution as the limiting 
distribution. If kn is unbounded, then the standard normal distribution is 
used for the appropriately standardized test statistic. 

The proposed test is simple and straightforward: it does not require esti- 
mation of /e(0) and the asymptotic null distribution of the test statistic is 
tractable. One might consider a Wald-type test as a natural choice for ^2 k„ 
under question. However, this would require estimation of /e(0) and would 
therefore be less desirable. 

For purpose of comparison, we consider the following likelihood-type-test. 
We use the objective function in (3.2) in place of the likelihood function and 
define the test statistic as 

{71 n 
Y^PriY, - ijx,) -Y^PriY, - if^k^Xi - i^,kT2,kATi,X,)) . 
i=l i=l J 



As in the parametric case, the likelihood-type test also requires estimation of 
/e(0) [23]. Following [1, 2], we simulate the null distribution via the residual 
bootstrap. 
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3.1. Heteroscedastic errors. Heteroscedastic error models are of great 
interest in quantile regression. The Rao-score-type test is applicable to a 
scale family of linear heteroscedastic models such that 

(3.4) Yi = e'^X, + siT„Xi)e,, 

where s(Tj, Xj) > is a scale function that is consistently estimable and the 
ei are i.i.d. 

A heuristic argument for this is as follows. The test only requires fitting 
the null model and the distribution of the test statistic relies only on the 
quantile estimator under the null. Therefore, any quantile estimator with 
the same first order representation under the null leads to the same lim- 
iting distribution for the test statistic. This implies that if we can obtain 
a -y/n-consistent estimator a, = s(Ti,Xi) and if the weighted quantile es- 
timator minimizing X^iLi ^i~^Pr(^i — d'^Xi) is first-order equivalent to the 
quantile estimator with s{Ti,Xi), then we can apply the Rao-score-type 
test to {{Ti, Xi)/ai,Yi/ai) as if the rescaled observations were i.i.d. A pri- 
mary example of interest is s(Tj,Xj) = {Ti,Xj-)j. We can obtain a ^/n- 
consistent estimator -y from the absolute residuals of median regression and 
Koenker and Zhao [25] showed that the weighted quantile estimator mini- 
mizing ^27=1 PriXi — Xi) with di = {Ti,Xj)^ is first-order equivalent 
to the quantile estimator with 7. For general s{Ti,Xi), Zhao [36] showed that 
s{Ti,Xi) can be consistently estimated using a nearest neighbor approach 
and that the resulting weighted quantile estimator is first-order equivalent. 
On the contrary, the likelihood-ratio-type test requires the alternative model 
to be fitted; hence, whether (or how) it can be applied to (3.4) is unclear. 

3.2. Test consistency. Both tests encounter a problem of inconsistency 
when kn is fixed since the constancy of the coefficient functions is tested via 
the constancy of the approximating polynomial splines. To be specific, we 
consider a set of functions L^,^ such that 

^kr,={q(t,x) GG\q{t,x) =^x - Rg^kAt^x)}, 

where ^j^^ = {iijO^)"^ and the Rg^knit^^) ^'I's given in Lemma A.l(ii). We 
note that L^^ is a subset of Q for which the approximations in the given 
spline space are functions of x with constant coefficients. When kn = ko, 
testing (3.1) means testing Ho:qT-{t,x) G L^y versus Hi:qr{t,x) G ^\LfcQ. 
Note that C C Lfc^, while Lfc(,\£ / 0. Also, any element of hk^XC is a null 
model and the proposed test will not be able to reject Hq, even as n— > 00. 
The set I^ko\^ constitutes a class of nonconstant functions in t associated 
with the fixed number of knots k^, against which the proposed test is not 
consistent. 

When kn is allowed to grow with n, we can consider a sequence {lLfc„} with 
L* = lim inf„^oo Lfc^ . With respect to this sequence, (3.1) can be restated 
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as Hq : qr{t, x) G L* versus Hi : qr{t, x) G Q\1L* . The following lemma shows 
that the test is consistent for all alternatives when A;„ increases. 

Lemma 1. Under Condition 1, L* = £. 

Fortunately, these issues do not pose a real concern in practice as the knots 
are chosen adaptively. For the null, the number of adaptively chosen knots 
will not grow with n and the results of [23] ensure that the proposed test 
is level-appropriate. For the alternative, the number of adaptively chosen 
knots will grow with n and we have the above lemma for consistency of the 
test along with Theorem 2 to ensure that the test is level-appropriate. 

4. Empirical results. 

4.1. Monte Carlo studies. We conducted Monte Carlo studies for the 
following i.i.d. and heteroscedastic error models in order to evaluate the 
proposed test relative to the likelihood-ratio-type test: 

Ml: y, = -ii+b{ti)xi + {ei-F~\T)), 

M2: = 72 + hxi + {-f^ti + 74Xi)(e, - ^"^(0.5)), 

where -Fe(") denotes the distribution function of the e^, the Xi are truncated 
standard normal variables, that is {xi = sign{zi) min(c, \zi\) with Zi standard 
normal and c > 0), the 7j are constants, ~ C/(0.5, 1.5) for all i and ti is 
independent of xi and for all i. Three different distributions were used 
for the e^: a standard normal, a with 1 degree of freedom and a t with 3 
degrees of freedom. Ml represents a homoscedastic error model that is a null 
model when h{ti) = h and an alternative model when h{ti) is a nonconstant 
function of t. Four different alternatives were considered for b{ti): a linear, 
quadratic, sine and log function of t. M2 represents a heteroscedastic error 
model that is a null model at the conditional median (r = 0.5). 

The Monte Carlo studies were based on 500 data sets, each a sample of 
size 200. We used piecewise linear splines with an adaptively chosen num- 
ber of uniform knots to fit a model to each simulated sample. The number 
of knots was chosen by the adapted SIC in (2.4). This closely resembles 
standard practice, in which the order of splines and knot placements are 
chosen nonadaptively, but the number of knots is chosen adaptively. Piece- 
wise linear splines were used because they provided reasonable fits for the 
various alternatives considered here. For the likelihood-ratio-type test, 200 
bootstrapped samples were used to simulate the null distribution of the test 
statistic. The performance of the tests was measured by power at the sig- 
nificance level a = 0.05. As the knots were chosen adaptively, powers under 
the null, that is the type I error rates, were expected to be higher than the 
nominal significance level. 
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Table 1 

Power at a = 0.05 when errors are i.i.d. (Ml). The parenthetical numbers in the second 
column indicate the densities of the respective error distributions at the rth quantile. 
"RS" and "LR" indicate the Rao-score- and likelihood-ratio-type tests 

Error constant linear quadratic sine log 

T distribution RS LR RS LR RS LR RS LR RS LR~ 

0.062 0.072 an 080 069 O80 084 O90 75 086 
0.040 0.062 1 1 1 1 1 1 1 1 
0.054 0.068 0.96 0.98 0.94 0.99 0.99 1 0.98 0.99 

0.070 0.130 0.39 0.64 0.4 0.64 0.43 0.75 0.46 0.65 
0.082 0.164 0.16 0.33 0.16 0.31 0.19 0.41 0.17 0.32 
0.056 0.130 0.40 0.62 0.49 0.70 0.50 0.78 0.55 0.72 



Table 1 summarizes the results for the homoscedastic error model. The 
likelihood-ratio-type test was slightly more powerful, but less likely to attain 
the nominal significance level, particularly when r = 0.9. The Rao-score- 
type test performed reliably with comparable power. Inherent to quantile 
regression, the performances of both tests depend on the error density at 
the quantile under consideration; power rises with increasing error density. 

Table 2 summarizes the results for the heteroscedastic error median null 
model. The Rao-score-type test attains type I error rates close to the nominal 
significance level. The likelihood-ratio-type test was administered as if the 
errors were i.i.d. since it was not clear how to simulate the null distribution 
of the test statistic via the bootstrap. Quite predictably, it did not perform 
well. 

4.2. Example. We illustrate the prescribed methodology using the forced 
expiratory volume (FEV) data from [29]. Measurements of FEV in liters, 
age (T) in years, height (if) in inches, sex (5 = 1 for boys/S = for girls) 
and smoking status {SM = 1 for a current smoker / SM = otherwise) were 
collected for 654 children aged 3-19 who participated in the Childhood Res- 
piratory Disease Study. 

Of interest are the effect of cigarette smoking on FEV and the lower con- 
ditional quantiles of FEV as a gauge for poor pulmonary function. The pro- 
posed methodology enables us to investigate both research questions while 
allowing the coefficients to vary with age to reflect the possibility of age- 
dependent covariate effects and the lower quantiles to be related to the 
covariates in ways different from those for the mean. We modeled the first 
quartile ((70.25) and the median (go.5) as follows: 

(4.1) Qr = Poit) + Piit) xS + p2it) xH + f33{t)x {H X S) + ^t) X SM. 

We tried both piecewise linear and quadratic splines with a set of eleven 
equispaced knots, using the stepwise knot selection algorithm described in 



0.5 Z (0.40) 
(0.67) 
^t, (0.64) 

0.9 Z (0.18) 
TfX? (0.09) 
4^*3 (0.18) 
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Table 2 

Type I Error Rate at a = 0.05 when errors are heteroscedastic 
(M2 with T = 0.5) 





Z 


— V? 


;^*3 


RS 


LR 


RS LR 


RS LR 


0.076 


0.378 


0.084 0.432 


0.056 0.324 



Intercept Height by Sex Interaction 




5 ID 15 5 10 15 

Fig. 1. Coefficient estimates by piecewise linear splines. 



Section 2.3. Both splines provided similar fits, while the fit with piecewise 
linear splines was preferred by the adapted SIC. Hence, we shall discuss the 
results based on piecewise linear splines. 

The estimates of /34(t) were small negative constants at both quantiles. 
This implies that the effect of smoking on FEV is not age dependent at these 
quantiles. This is not surprising since the current smokers in the data set 
were close to each other in age. 

/3o(i)i f^2(t) and (33{t) were estimated as varying with age (see Figure 1). 
At both quantiles, the estimated intercept increases linearly, implying that 
FEV increases linearly with age. The height coefficient for girls was esti- 
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mated as a constant for qo_25, while knots were chosen at 5 and 11 for (70.5- 
For the coefficient of height by sex interaction, knots were chosen at 6 for 
Qo.25 and at 8 for go. 5- The height coefficient estimates for boys in the bottom 
right panel were obtained as P2{t) + Psit)- 

We focus on the estimates for ages from 5 to 15 in Figure 1 since little data 
is available outside the age range. The nonparallel intercept estimates and 
nonoverlapping coefficient estimates imply a different regression relationship 
at the first quartile than at the median. More specifically, the estimated 
intercept does not grow as rapidly at the first quartile as at the median. 
A similar pattern is shown in the height coefficient estimates for boys. This 
suggests that in boys, FEV does not increase with age as rapidly in the lower 
tail of its conditional distribution as in the middle, implying that in boys, 
the difference between height-adjusted typical and poor pulmonary function 
increases with age. It manifests through increasing variability in FEV over 
that age range, while variability in height does not necessarily increase (see 
Figure 2). This structural heterogeneity would not have been revealed by an 
ordinary regression. 
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We also tested for age-dependent effects of the covariates using the pro- 
posed Rao score test. We used hnear sphnes with the same knots selected 
for the regressions. We rejected the constancy hypothesis with p- values near 
zero for both quantiles. 

While FEV is known to increase with age, it is less clear which other 
coefficients change with age. In particular, do the coefficients of height and 
height by sex interaction change? To address this question, we detrended 
the data by subtracting age-and-gender-specific sample quartiles from the 
observations of respective age-gender groups. We then modeled the condi- 
tional quartiles of the detrended data as a function of height and height 
by sex interaction without an intercept. We used piecewise splines with the 
same knots selected for the regressions. We rejected the constancy hypoth- 
esis with a p- value < 0.006 for (70.25 and a p- value near zero for go. 5- The 
results agreed with the nonconstant coefficient estimates obtained from the 
regressions. 

APPENDIX 

For simplicity, proofs are provided for uniform knots and p = 2, the sub- 
script kn being suppressed if possible. The same proofs will work for quasi- 
uniform knots and/or p> 2. We adopt the following notation: 

Gn = (vr(ri), . . . , 7r(r„,))(vr(r,), . . . , 7r{Tn)f , 

-^n = (n(T'j, Xj), . . . , n(r„, x„))p^^ xn, h„ = xj^x^, 
G = E(7r(r)7r(r)T), H = E(n(r,x)n(r,x)T), ^.. = ^^.,jro, 

Zi = U-^/^U{Ti,Xi), 

1/2 

where denotes the Moore inverse of Hn . Also, for a matrix D, 

let D(j j) denote its element at the ^th row and jth column. For example, 

^n{aN+c,bN+d) = Yl,i'=l Bc{Ti)Bd{Ti)Xi(^a)Xi{b) and il{aN+c,bN+d) = 

E{Bc{T)Bd{T)X(^^)X^^)) for o, 6 = 0, 1, 2 and 1 < c, d < /c„ + m + 1. 

We first present lemmas that are necessary to prove the theorems. Lemma 
A.l follows directly from Corollary 6.21 of [30] and the proofs of Lem- 
mas A.2-A.4 follow lines of argument similar to those used in the proofs 
of lemmas in [5]. Lemma A. 5 follows directly from Lemma A. 4 and Condi- 
tion 3. Therefore, we omit their proofs. 

Lemma A.l. Assume Conditions 1 and 3. For some constants Wi and 
Wi that depend on only (hut not all of) m, Wq, p and M : 
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(i) there exist6jTTk„{t), j = 0,...,p, such that j3j{t) = OjTTk„{t) -rj^k^it) 
and supig[o^i] \rj,k,M)\ < Wik'"' ; 

(ii) there exists 9 '^^^Ilk„{t,x) such that qT-{t,x) = 6'i^^Ilf^^{t,x) — Rg k^{t,x) 
and sup(j_^.)g[o,i]xRP+i \Rg,k„it,^)\ < ^i^n''- 

Lemma A. 2. Assume Conditions 1-3. For all kn, the eigenvalues of 
knG and /c^H are bounded and depend only on m. 

Lemma A. 3. A SSUTflC CoTlditiOTlS 1^3. Suppose lilllj^^oo knfi^ ^ = for 
some < 6 <1. Except on an event which depends on {{Ti,Xi), {T2,X2), ■ ■ ■ , 
{Tn,Xn)} and whose probability tends to zero with increasing n, \iln{aN+c,bN+d)/ 

n — 'H.(^aN+c,bN+d)\ = Op{cnn~^^'^ kn^"^) , Uniformly in a,b,c and d, where {c„} 
is a nondecreasing sequence of positive numbers for which c^/nk~^ and 
cl/log kn^oo. 

Lemma A. 4. A ssume Conditions 1—3. Suppose lini^_>oo knu^-^ = for 
some < (5 < 1. Then the eigenvalues of knGn/n and fc„H.„/?i are bounded 
in probability. 

Lemma A. 5. Assume Conditions 1-3. Suppose liuin^oo knn^~^ = for 
some < 6 < 1. Then with probability one, maxi<j<„ \Zi\ = 0{\/kn/n). 

Lemma A. 6. Under the conditions of Theorem 1, {6 — 0)'^H„(^ — 0) = 

Op{kn) ■ 

Proof. Since {§ - 0)'^H„(0 -6) = {6* - e*)^{e* - §*), it suffices to 
show that \0 - p = Op{kn). As 6 = arg min0gigPfe„ Yd=iPr{<^i - {0 - 
e y^Zi - Ri) from (2.2), with probabihty 1 we have 

n n 

(A.l) y^p^{e,-{e* -e*fZi-Ri)= inf y^pr{ei-e^Z,-Ri). 

Let Ql = {0\6 G W^n , \0\ < Lk}/^} for some constant L > 0. Then the proof 
is reduced to showing that a solution to the optimization problem in (A.l) 
lies in for a sufficiently large L. The following lemma and the convexity 
of the minimization problem (A.l) complete the proof. □ 

Lemma A. 7. Assume the conditions of Theorem 1. 

(i) For any sequence {Ln} satisfying I < Ln < kn°^^^ for some < ?yo < 
(r-l/2)/(2r + l), we have 



sup /c„ 



J2[pr{e^ - Zfe - Ri) - - R^) + Zje{T - I{ei < 0)) 



1=1 
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E{pr{ei - Zje - Ri) - prie, - R,))] 



opil). 



(ii) For any e > 0, there exists L := (sufficiently large) such that as 
n — > oo, 

p\k-' ( inf J2i^{pr{ei - Zje - Ri) - pria - R,))] 



,\e\=Lk'Ju=i 



J2Z^iT-I{e^<0)) 



1=1 



>U >l-e. 



The proof of Lemma A. 7 uses the results of Lemmas A.4-A.5 and a par- 
tition of the parameter space It apphes Bernstein's inequahty to bound 
the probabihties over the partition. The method is similar to that used in 
Lemmas 3.2 and 3.3 of [14] and so we omit the details. 

By Lemma A. 7, there exists some L := for any e > such that as 
n — > oo, we have 

{n n ^ 

inf Y^Pri^i - Zje - R,) >Y,Pr{e^ - R^) \ > I - e. 
\e\=Lki^^ i=i i=i J 

By Corollary 25 of [8], we have that as n — > oo, 

{n n \ 

inf 5]p.(ei - Zfe -R^)>Y. pM -R,)\>l-e. 
\e\>Lkl'\=i i=i J 

By (A.2), (A.l) implies that - 6**1 < LkH'^] > 1 - e. 

Proof of Theorem 1. First, it follows from Lemma A.l(i) that 

1 n r) n r) n 



n 



i=l 



n 



i=l 



n 



i=l 



n 



n 



As kn ~ ?7,^/(^''+^\ by Lemma A. 4, it suffices to show that \6j — 9j\'^ = 
Op{k'^/n). By Lemmas A. 4 and A. 6, \0 — 0p = Op{k'^/n), which implies 
Op{kl/n). □ 



that lOj - 6'j|2 



Proof of Lemma 1. As £ C L* by definition, it suffices to show that 
L* C C. Suppose that L*/£ / 0. Let q*{t,x) denote an element of L*/£ 
and let {kn*} denote a subsequence of {kn} such that q*{t,x) S L^^^ for all 
kn* - Consider Rq*^k„, {t, x), where Rq*^kn» (*) is given in Lemma A.l(ii) with 
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respect to q*{t,x). We note that q* G L*/£ implies limk^^^oo Rg,kn,{'t,x) / 0. 
This contradicts Lemma A.l(ii). □ 

Proof of Theorem 2. Let Q = {Q^n'^^y^'^T2{Ti,Xi). Define = 1/a 
lA=iCi^r{ei) and bk„ = l/crELi Oiv^rCci - (|i -ii)^Xi) - (friei))}. Then 
Wk„ = ak„ + bk„. We have l°fcnl'-(Pfcn-P-i) — ,D ^^q^ ^^^^ Theorem 4.1 

of [28] by replacing the continuity condition on (pr{u) in the theorem with 
Condition 2. As \aj^bkj < \akj\bkj and Pfc„ - p - 1 = (p + 1)(A:„ + ?n), it 
suffices to show that \bk„ \ = Op{\/l{^). The following lemma completes the 
proof. □ 

Lemma A. 8. Assuming the conditions of Theorem 2, \bk^ \ = Op{\fk^). 

Proof. First, Fe{t) - Fe{0) = fe{0)t + o{\t\) and E{ipr{s + t) - ipr{s) f = 
/e(0)i + o{\t\) for sufficiently small Also, from Lemma A. 5, maxj|^j| = 
0{\/kn/n) with probability one. 

Consider a function 6(r/) = l/cX^iLi CilVrlej — xfrj/^/n) — ipri^i)} for 
some rj G ]R^^+^. For b(^j-^{r]), j = 1, . . . ,{p + l){kn + m), we have 

n 
i 

and l/fjmaxj |Ci{j){'/'r(ei — xjr]/ y/n) — ipr{^i)} \ = 0{^kn/n). It follows from 
Bernstein's inequality that for all L > 0, e > and any |r/| < L, 

P{|6(j)(r7) - > e} < 2exp(- V^£V(2A:n + 2v^/3)). 

As the number of 6(j)(?7)'s to consider is (p + l){kn + m) and by the arbi- 
trariness of L > and e > 0, we have 

P{ sup \b{5)- E{b{r}))\ >e} 

< 2exp(-^eV(2fcn + 2yfc;;/3) + log((p + l)(fe„ + m))). 

Under the conditions of Theorem 2, nk~'^ / (log /c„)^ — > oo and the right-hand 
side converges to zero as n ^ oo. On the other hand, \E{b{if]))\ = 0{^/k^Jn), 
except for an event which depends on { (Ti, Xi), (r2, X2), . . . , (T^, X„) } and 
whose probability tends to zero with increasing n. Since |0ifc„| = Op{l/^/n)^ 
the lemma follows. □ 
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